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Communication shapes sensory response in multicellular networks 
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Collective sensing by interacting cells is observed in a variety of biological systems, and yet a 
quantitative understanding of how sensory information is collectively encoded is lacking. Here we 
investigate the ATP-induced calcium dynamics of monolayers of fibroblast cells that communicate 
via gap junctions. Combining experiments and stochastic modeling, we find that increasing the 
ATP stimulus increases the propensity for calcium oscillations despite large cell-to-cell variability. 
The model further predicts that the oscillation propensity increases not only with the stimulus, but 
also with the cell density due to increased communication. Experiments confirm this prediction, 
showing that cell density modulates the collective sensory response. We further implicate cell-cell 
communication by coculturing the fibroblasts with cancer cells, which we show act as “defects” 
in the communication network, thereby reducing the oscillation propensity. These results suggest 
that multicellular networks sit at a point in parameter space where cell-cell communication has a 
significant effect on the sensory response, allowing cells to simultaneously respond to a sensory input 
and to the presence of neighbors. 


I. SIGNIFICANCE STATEMENT 

Cells routinely sense and respond to their environment, 
and they also communicate with each other. Communi¬ 
cation is therefore thought to play a crucial role in sens¬ 
ing, but exactly how this occurs remains poorly under¬ 
stood. We study a population of fibroblast cells that 
responds to a chemical stimulus (ATP) and communi¬ 
cates by molecule exchange. Combining experiments and 
mathematical modeling, we find that cells exhibit cal¬ 
cium oscillations not only in response to the ATP stim¬ 
ulus, but also in response to an increase in cell density. 
To confirm that the density dependence is a result of in¬ 
creased cell-cell communication, we combine the fibrob¬ 
lasts with cancer cells, which we show have weakened 
communication properties. The oscillations indeed de¬ 
crease with the fraction of cancer cells. Our results show 
that when cells are together, their sensory responses re¬ 
flect not just the stimulus level, but also the degree of 
communication within the population. 


II. INTRODUCTION 

Decoding the cellular response to environmental per¬ 
turbations, such as chemosensing, photosensing, and 
mechanosensing has been of central importance in our 
understanding of living systems. Up to date, most studies 
of cellular sensation and response have focused on single 
isolated cells or population averages. An emerging pic¬ 
ture from these studies is the set of design principles gov¬ 
erning cellular signaling pathways: these pathways are 
organized into an intertwined, often redundant network. 
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and the network architecture is closely related with the 
robustness of cellular information processing [U [2] . On 
the other hand, many examples suggest that collective 
sensing by many interacting cells may provide another 
dimension for the cells to process environmental cues [3] . 
Quorum sensing in bacterial colonies [4] and cyclic adeno¬ 
sine monophosphate (cAMP) signaling in Dictyostelium 
discoideum populations [5] generate dramatic group be¬ 
haviors in response to environmental stimuli. Various 
types of collective sensing have also been observed for 
higher level multicellular organisms, such as in the olfac¬ 
tion of insects [ 6 ] and mammals [ 7 ] , the glucose response 
in the pancreatic islet [8], and the visual processing of 
animal retinal ganglion cells mm- These examples sug¬ 
gest a fundamental need to revisit cellular information 
processing in the context of multicellular sensation and 
responses, since even weak cell-to-cell interaction may 
have strong impacts on the states of multicellular net¬ 
work dynamics [9]. In particular, we seek to examine 
how the sensory response of cells in a population differs 
from that of isolated cells, and whether we can tune be¬ 
tween these two extremes by controlling the degree of 
cell-cell communication. 

Previously, we have described the spatial-temporal dy¬ 
namics of collective chemosensing of a mammalian cell 
model system muig. In this model system, high den¬ 
sity mouse fibroblast cells (NIH 3T3) form a monolayer 
that allows nearest-neighbor communications through 
gap junctions m- When extracellular ATP is deliv¬ 
ered to the cell monolayer, P2 receptors on the cell mem¬ 
brane trigger release of the second messenger IP3 (Inos¬ 
itol 1,4,5-trisphosphate) [14]. IP3 molecules, when bind¬ 
ing to their receptors at the membrane of endoplasmic 
reticulum (ER), open ion channels to allow influx of Ca^+ 
from the ER to the cytoplasm. The fast increase in cy¬ 
tosolic calcium concentration is followed by a slow de¬ 
crease when Ca^+ are pumped back to the ER [15] . This 
simple picture, however, is complicated by the nonlinear 
feedback between Ca^+ and ion channel opening prob- 
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ability, which leads to rich dynamic behaviors such as 
cytosolic calcium oscillations m- In the situation of 
collective ATP sensing, we have found that gap junc¬ 
tion communications dominate intercellular interactions 
m- Furthermore, these short-range interactions propa¬ 
gate and turn the cell monolayer into a percolating net¬ 
work [T2j. These characteristics make the system ideal 
for studying how sensory responses are modulated by ex¬ 
tensive communication in multicellular networks. 

Here we use this model system to examine how cell-cell 
communication affects collective chemosensing. Combin¬ 
ing experiments with stochastic modeling, we find that 
cells robustly encode the ATP stimulus strength in terms 
of their propensity for calcium oscillations, despite signif¬ 
icant cell-to-cell variability. The model further predicts 
that the oscillation propensity depends not only on the 
stimulus, but also on the density of cells, and that denser 
monolayers have narrower distributions of oscillation fre¬ 
quencies. We confirm both predictions experimentally. 
To verify that the mechanism behind the density depen¬ 
dence is the modulation of cell-cell communication, we 
introduce cancer cells (MDA-MB-231) into the fibrob¬ 
last cell monolayer. As we show, MDA-MB-231 cells act 
as “defects” in the multicellular network, as they have 
distinct calcium dynamics when compared with the fi¬ 
broblasts due to reduced gap junction communication 
[HHIl]. We find that the oscillation propensity of the 
fibroblasts decreases as the fraction of cancer cells in¬ 
creases, confirming that the sensory response is directly 
affected by the cell-cell communication. Our findings in¬ 
dicate that cells’ collective response to a sensory stimu¬ 
lus is distinct from their individual responses, and that 
in this case the response simultaneously encodes both 
the stimulus strength and the degree of communication 
within the population. 


III. RESULTS 

In order to study the sensory responses of a multicel¬ 
lular network, we use single channel microfluidic devices 
and deliver ATP solutions to monolayers of fibroblast 
(NIH 3T3) cells (see Appendix A). The stimuli solutions 
of ATP concentrations varying from 0 jaM to 200 jaM 
are pumped though the channel at rate of 50 /il/min, 
ensuring minimal flow perturbation to the cells and fast 
delivery of ATP across the field of view. During this 
time we record the calcium dynamics of each individ¬ 
ual cell loaded with fluorescent calcium indicator (Fluo- 
4, Life Technologies) at 4 frames/sec (Hamamatsu Flash 
2.8, Hamamatsu). Our analysis uses 400 seconds of data 
from each recording selected such that ATP arrival is at 
approximately 50 seconds. 

We modulate the degree of communication in two 
ways. First, we vary the cell density. Smaller cell densi¬ 
ties correspond to larger cell-to-cell distances, which we 
expect to reduce the probability of forming gap junctions. 
Second, we coculture the fibroblasts with breast cancer 


(MDA-MB-231) cells in the flow channel (see Appendix 
A). As we later show, MDA-MB-231 cells have reduced 
communication properties and therefore act as defects in 
the multicellular network. To distinguish the two cell 
types, MDA-MB-231 cells are pre-labeled with red flu¬ 
orescent dye (Cell Tracker Red CMTPX, Life Technolo¬ 
gies) . Varying cell density and the fraction of cancer cells 
allows us to control the architecture of the multicellular 
network over a wide range. 

Figureshows the composite image of a high-density 
cell monolayer with cocultured fibroblast and cancer 
cells. In this example, MDA-MB-231 cells make up a 
fraction Fc = 15% of the total population which has a 
total cell density of pT = 2500 cells/mm^. At this den¬ 
sity, each cell has an average of 6 nearest neighbors from 
which extensive gap junction communication is expected. 
After identifying cell centers from the composite image 
(see Appendix A), we compute the time-dependent av¬ 
erage fluorescent intensity near the center of each cell 
which represents the instantaneous intracellular calcium 
concentration at the single cell level. 


A. Collective response to ATP stimuli 

Typical responses of cells to excitation by ATP are 
shown in Fig. We see that, on average, higher con¬ 
centrations of ATP trigger larger increases in calcium 
levels. Cell-to-cell variations are significant; for example, 
response times as well as subsequent calcium dynamics 
of individual cells vary dramatically amongst cells in the 
same multicellular network. In many cells, the initial 
calcium increase is followed by transient calcium oscil¬ 
lations. We quantify the oscillation propensity by com¬ 
puting the fraction of non-oscillating cells Fn using a 
peak-finding algorithm (see Appendix B). We see in Fig. 
EP that higher concentrations of ATP cause a larger per¬ 
centage of cells to oscillate, and thus a smaller Fat. 

The period of the oscillation is characterized by the 
inter-spike interval (ISI), which has been proposed to dy¬ 
namically encode information about the stimuli [20l |2T] . 
To investigate the characteristics of ISI in the context col¬ 
lective chemosensing, we systematically study the statis¬ 
tics of the ISI from 30,000 cells. Figure shows the 
histogram (event counts) of ISI values, normalized by 
the number of cells, of a typical experiment where the 
ATP concentration is 50 pM. We see that the distribu¬ 
tion is broad, which underscores the high degree of cell- 
to-cell variability in the responses. Figure HP summa¬ 
rizes the distribution at each ATP concentration using a 
box-and-whisker plot. We see that there is no significant 
dependence of the ISI on the ATP concentration. This is 
at odds with a familiar property of calcium oscillations, 
termed frequency encoding, in which the oscillation fre¬ 
quency (or ISI) depends on the strength of the stimulus 
pTl [2Ql [22l [23]. However, we will see in the next section 
that the lack of a dependence here is likely due to the 
high degree of cell-to-cell variability. 
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FIG. 1: Calcium dynamics of cell monolayer in response to 
extracellular ATP. (A) Composite image showing the multi¬ 
cellular network of cocultured fibroblast (NIH 3T3) and breast 
cancer cells (MDA-MB-231). Red: MDA-MB-231. Green: 
fluorescent calcium signal for all cells (MDA-MB-231 and NIH 
3T3). (B) Normalized fluorescence intensity profiles of one 
typical experiment for each ATP concentration tested. Blue: 
randomly selected single cell calcium responses. Red: average 
intensity profiles of all cells in each experiment. All time series 
begin approximately 50 seconds before arrival of ATP stim¬ 
uli. Intensity profiles of individual cells have been rescaled to 
[—1, 1]. (C) Fraction of non-oscillating cells Fn as a function 
of ATP concentration at fixed cell density {pr = 1200 ± 200 
cells/mm^) and cancer cell fraction (15% ± 6%). Error bars: 
SEM for n > 4. (D) ISI event counts normalized by number 
of cells for only NIH 3T3 cells. (E) Average experimentally- 
measured ISI values of NIH 3T3 cells at varying ATP concen¬ 
trations at fixed cell density {pr = 1200 ±200 cells/mm^) and 
cancer fraction {Fc = 15% ± 6%). (E) ISI cross-correlation 
as a function of topological distance. Data from experiments 
with 50 pM ATP, at fixed cell density {pr = 1400 ± 400 
cells/mm^) and cancer fraction (Fc = 20% ±5%). Error bars 
show standard deviation from five experiments. 


Finally, we characterize the spatial correlations of 
the ISI within the monolayer by computing the cross¬ 
correlation function Cisi{d) as a function of topologi¬ 
cal distance d between cells. In particular, for each ex¬ 
periment, we first identify all oscillatory cells and com¬ 
pute the average ISI Ti for each cell i. We then define 
6Ti =Ti- (Ti), and Cisi{d) = {STiSTj)D,,=d/{ST^), 
where Di^s the topological distance between cells i and 
j. Figure]^ shows that Cisi{d) falls off immediately for 
d > 0. This is surprising, since the cells experience iden¬ 
tical ATP stimuli, and one might hypothesize that com¬ 
munication between cells would result in the ISI values 
for nearby cells being correlated. However, as described 
next, evidence from mathematical modeling suggests that 
this correlation is removed by the cell-to-cell variability. 


B. Stochastic modeling of the collective response 

To obtain a mechanistic understanding of the experi¬ 
mental observations, we turn to mathematical modeling. 
We develop a stochastic model of collective calcium sig¬ 
naling based on the work of Tang and Othmer [20l [24] . 
Their model captures the ATP-induced release of IP3, 
the IP3-triggered opening of calcium channels in the ER, 
and the nonlinear dependence of the opening probability 
on the calcium concentration, as schematically illustrated 
in Fig. The model neglects more complex features of 
calcium signaling observed in some cell types, such as up¬ 
stream IP3 oscillations [25l [26] and spatial correlations 
among channels [271128]. The model predicts that at a 
critical ATP concentration, the calcium dynamics transi¬ 
tion from non-oscillating to oscillating. However, it was 
previously only analyzed deterministically for a single cell 
[20l[24|. Therefore we extend it to include both intrinsic 
noise and cell-cell communication via calcium exchange 
(see Appendix C). We also explicitly include the dynam¬ 
ics of IP3, which has a constant degradation rate and a 
production rate a that we take as proportional to the 
ATP concentration. We simulate the dynamics using the 
Gillespie algorithm [29], and we vary the density pT of 
cells on a square grid, which modulates the degree of 
communication. 

Figure shows the dependence of Fn on a, where 
a is the model analog of the experimentally controlled 
ATP concentration. Consistent with the experimental 
findings in Fig. [^, we see that Fn decreases with a. In 
the model, the decrease is due to the fact that intrinsic 
noise broadens the transition from the non-oscillating to 
the oscillating regime. Thus, instead of a sharp transition 
from Fn = 1 to Fn = 0 as predicted deterministically, 
the transition occurs gradually over the range of a shown 
in Fig. 1^. Figure shows the dependence of the ISI 
on a in the model (see the green box plots). We see 
that the ISI decreases with a, which is expected since 
frequency encoding is a component of the Tang-Othmer 
model [20l[24]. Yet, this property is not consistent with 
the experimental observation in Fig. [^, where the ISI 
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shows no clear dependence on ATP concentration. Fur¬ 
ther, Fig. shows the dependence of the correlation 
function Cjsi on the topological distance d in the model 
(see the green curve). We see that Cjsi decreases gradu¬ 
ally with d, indicating nonzero spatial correlations in the 
ISL This feature is again inconsistent with the experi¬ 
mental findings seen in Fig. ET- 

Motivated by the high level of cell-to-cell variability ev¬ 
ident in Figs. and [IP> we hypothesize that cell-to-cell 
variability is responsible for these discrepancies between 
the model and the experiments. Indeed, inspecting the 
ISI histogram from the model reveals a very narrow dis¬ 
tribution of ISI values, as seen in Fig. (green bars), 
which is in contrast to the broad distribution observed ex¬ 
perimentally in Fig. HP- To incorporate cell-to-cell vari¬ 
ability, we allow the model parameters to vary from cell 
to cell. Lacking information about the susceptibility of 
particular parameters to variation, we allow all model pa¬ 
rameters to vary by the same maximum fold change (i.e., 
parameters are sampled uniform randomly in log space). 
The maximum fold change F is found by equating the 
variance of the resulting ISI distribution with that from 
the experiments, which yields the value F = 2 (see Ap¬ 
pendix C). As seen in Fig. (blue bars), the resulting 
ISI distribution is consistent with that observed in Fig. 
[2p, both in width and in shape. 

We see in Fig. (blue box plots) that including cell- 
to-cell variability in the model severely weakens the de¬ 
crease of the ISI with d, therefore agreeing with the ex¬ 
perimental results shown in Fig. HP (with variability, a 
is defined as the mean of the a values sampled for each 
cell). We also see in Fig.(blue curve) that variability 
removes the correlations Cjsi for d > 0, which is con¬ 
sistent with the immediate fall-off observed experimen¬ 
tally in Fig. HP- Importantly, even with variability, the 
decrease of F^ with a seen in Fig.[^ persists, as demon¬ 
strated in Fig. [^. This decrease remains consistent with 
the experimental observation in Fig. HP- Indeed, variabil¬ 
ity significantly broadens the range of a values over which 
the transition occurs, as expected (compare Fig. to 
Fig. 1^), which is consistent with the broad range over 
which the transition occurs experimentally (Fig. EP)- 


C. Effects of communication on the sensory 
response 

Having validated the model, we now use it to make 
predictions about the effect of cell-cell communication 
on collective calcium dynamics. Communication in the 
model is controlled by cell density, with higher density 
leading to more cell-to-cell contacts and thus a higher 
degree of communication. Therefore we first investigate 
the dependence of the oscillation propensity on the cell 
density. Fig. IP shows F/v as a function of both cell den¬ 
sity pT and the ATP-induced IPS production rate a. We 
see that the fraction of non-oscillating cells transitions 
from Fn = 1 to F/v = 0 as a function of a and that there 
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FIG. 2: Model development and validation. (A) Schematic 
of the model. ATP stimulates IPS release at rate a, and 
IPS acts jointly with Ca^^ to open calcium channels (positive 
feedback), while further Ca^^ binding closes channels (neg¬ 
ative feedback). Communication is modeled via diffusion of 
Ca^^ between adjacent cells. (B) Fraction of non-oscillating 
cells Fat as a function of ATP-induced IPS production rate 
a. (C) Interspike interval (ISI) decreases with a (green). The 
decrease is severely weakened by cell-to-cell variability (blue). 

(D) ISI cross-correlation as a function of topological distance d 
(green). Cell-to-cell variability removes correlations for d > 0. 

(E) Distribution of ISI values (green). Cell-to-cell variability 
significantly broadens distribution (blue). (F) Fn versus a 
with cell-to-cell variability. In B and F, error bars are SEM 
for n = 5 subsamples. 


is also a dependence of F^ on pT- Not only has intrinsic 
noise and cell-to-cell variability broadened the Fn tran¬ 
sition as a function of d, but evidently cell-cell communi¬ 
cation has introduced a dependence of the transition on 
cell density. The result is that at low d, Fn is everywhere 
large and is independent of cell density (Fig.|^ and left 
box in Fig. |^). However, at intermediate a, Fn is a de¬ 
creasing function of cell density (Fig.[^ and right box in 
Fig.[^). In this regime, increasing the degree of commu¬ 
nication causes more cells to exhibit oscillatory calcium 
dynamics (thus decreasing F/v), even with a fixed sensory 
stimulus a. At large a (beyond the range of Fig.|^), we 
have checked that the non-oscillating fraction is driven 
to low values as expected, and the density dependence of 
Fn diminishes. 

The prediction in Fig. is striking because it im¬ 
plies that cell-cell communication causes more cells to 
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FIG. 3: Model predictions. (A) Fraction of non-oscillating 
cells FV as a function of cell density pr and ATP-induced 
IPS production rate a. Left and right boxes correspond to B 
and C, respectively. (B) At small a, Fn is large and density- 
independent. (C) At intermediate a, Fn decreases with den¬ 
sity. In B-D, error bars are SEM for n = 5 subsamples. (D) 
Entropy of ISI distribution Hjsi increases with Fn- 


oscillate, even while cell-to-cell variability causes their 
ISI values to be spatially uncorrelated (Fig. [^). There¬ 
fore, we wondered whether communication would have 
an effect on the width of the ISI distribution in this 
regime. The width, or more generally the amount of 
uncertainty contained in the ISI distribution, is char¬ 
acterized by the entropy. For a continuous variable x, 
the entropy becomes the differential entropy, defined as 
Hjsi = — f p{x) log p{x)dx, where p{x) is the probabil¬ 
ity density. As seen in Fig. [^, the entropy of the ISI 
distribution increases with Fn- This indicates that as 
communication decreases the non-oscillating fraction, it 
also narrows the distribution of ISI values. 

We now test these predictions in our experimental 
system. To test our predictions about how the non¬ 
oscillating fraction of cells should depend on cell den¬ 
sity, we measure Fn as a function of pT for various ATP 
concentrations. We see in Fig. UK that with no ATP, 
Fn is large at both low and high densities, and there 
is no statistically significant correlation between Fn and 
pT- Then, we see in Fig. that at intermediate ATP 
concentrations (10—100 yuM), Fn significantly decreases 
with pT (see Appendix B for pairwise statistical com¬ 
parison of particular density values). Finally, we see in 
Fig.|4p that at large ATP concentration (200 /iM), Fn 
is small at both low and high densities, and again there 
is no statistically significant correlation between Fn and 
pT- These results confirm the predictions in Fig. 

To test the prediction that the entropy of the ISI distri¬ 
bution increases with the non-oscillating fraction of cells, 
we measure Hjsi as function of Fn- As seen in Fig. iP, 
Hi SI increases with Fat, consistent with the prediction 
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FIG. 4: Experimental tests of model predictions. (A) Frac¬ 
tion of non-oscillating NIH 3T3 cells Fat as a function of cell 
density pr when stimulated by 0 pM ATP. Error bars: SEM 
for n > 4. ns, not significant. (B) As in A, but with inter¬ 
mediate concentrations 10—100 pM ATP. (G) As in A, but 
with 200 pM ATP. In A-G, the cancer cell fraction is hxed 
at Fc = 15% ± 6%. (D) Fn is positively correlated with the 
differential entropy of inter-spike intervals Hi si- Error bars 
represent standard deviation of 1000 bootstrap resampled re¬ 
sults, see Appendix D for more details. 


in Fig. Ip. This implies that cell density regulates the 
spectrum of the ISI response. In particular, it suggests 
that increasing the degree of cell-cell communication nar¬ 
rows the distribution of ISI, making the ISI values less 
variable across the cell population. We have also checked 
that the entropy of the distribution of cross-correlation 
values for nearest neighbors’ entire calcium trajectories 
Cnn dHHl] decreases as a function of cell density (see 
Appendix D). Cnn is not only a statistical characteri¬ 
zation of the collective cellular dynamics, but also may 
provide robust channel for cellular information encoding 
ElEolEI]. Together, these results imply that cell-cell 
communication has a significant effect on the collective 
sensory response. This finding is especially striking given 
the strong effects of cell-to-cell variability (Fig. HP and 
F). We conclude that the effects of communication ob¬ 
served here persist in spite of extensive cell-to-cell vari¬ 
ability. 


D. Effect of cancer cell “defects” 

We have seen that increasing cell density increases the 
propensity of cells to oscillate in response to an ATP 
stimulus. This behavior is consistent with our model. 
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which predicts that the mechanism is through increased 
cell-cell communication. However, it could be in the 
experiments that increasing the cell density introduces 
other effects beyond increased gap-junction communica¬ 
tion, such as mechanical coupling between cells or cou¬ 
pling to the substrate [32]. To modulate the communi¬ 
cation directly, we vary the fraction Fc of cancer cells 
with which the fibroblasts are cocultured, while keeping 
the density of all cells fixed. Since cancer cells are known 
to have reduced gap junction communication [TTHIS] , we 
expect the fraction of non-oscillating cells F/v to have the 
opposite dependence on Fc that it does on cell density 
(Fig.|^). 

First we investigate whether MDA-MB-231 cells in¬ 
deed have reduced communication in our system. Fig¬ 
ure shows several examples of single-cell calcium dy¬ 
namics for NIH 3T3 and MDA-MB-231 cells in a typical 
experiment. We see that both cell types exhibit imme¬ 
diate increases in cytosolic calcium levels at the arrival 
of ATP, but cancer cells typically show long relaxation 
times while fibroblast cells tend more often to exhibit 
oscillations after stimulation. These qualitative features 
are maintained across all ATP concentrations. Figure 
shows a comparison of the intercellular diffusion co¬ 
efficients in the two cell types, obtained from a fluores¬ 
cence recovery after photobleaching (FRAP) analysis [33| 
(see Appendix A). We see in Fig. that gap junction- 
mediated diffusion between MDA-MB-231 cells is signif¬ 
icantly weaker than between NIH 3T3 cells, consistent 
with previous reports niHis]. It is therefore evident that 
MDA-MB-231 cells can be treated as communication de¬ 
fects in the co-cultured multicellular network. Indeed, 
Fig. shows the spatial distribution of these defects in 
the monolayer. In Fig. [^, the mean ISI for each cell is 
shown in color, with non-oscillating cells in black. We see 
that cancer cells, labeled by white circles, are more likely 
to be non-oscillating, which is consistent with the quali¬ 
tative characteristics shown in Fig. IK- We have further 
quantified the distinction between the two cells types in 
Appendix B, where we show using the distributions of 
ISI values that oscillatory events are at least five times 
less likely to occur for the MDA-MB-231 cells. 

Having established that the presence of cancer cells re¬ 
duces the degree of cell-cell communication in the mono- 
layer, we now vary the fraction of cancer cells and mea¬ 
sure the oscillation propensity of the remaining fibrob¬ 
lasts. Fig. [^ shows the non-oscillating fraction of fi¬ 
broblasts Fn (blue bars) as a function of the cancer cell 
fraction Fc for a typical experiment at fixed cell density 
{pT = 1200 ± 200 cells/mm^). We see that Fn signif¬ 
icantly increases with Fc (see Appendix B for pairwise 
statistical comparison of particular Fc values). We also 
see that F/v for all cells (both fibroblasts and cancer cells, 
red bars) significantly increases with F^, and that, as 
expected, F/v is larger for all cells than for just fibrob¬ 
lasts. These findings imply that reduced cell-cell commu¬ 
nication decreases the propensity for calcium oscillations, 
which is consistent with the effects of varying cell density 



FIG. 5: Effects of cancer cell “defects” on collective response. 
(A) Typical fluorescence intensity profiles showing the cal¬ 
cium dynamics on the single cell level, where basal level in¬ 
tensity has been subtracted. For each cell, basal level inten¬ 
sity is estimated by averaging 100 seconds of its fluorescent 
intensity before ATP arrival. (ATP concentration = 50 pM, 
Pt = 2400 cells/mm^, Fc = 12%). (B) Fluorescence recov¬ 
ery after photo bleaching (FRAP) experiments confirm that 
MDA-MB-231 cells have weaker gap junction communication 
compared with NIH 3T3 cells (error bars: SEM for n > 100). 
See Appendix A for more details. (C) Spatial map of av¬ 
erage ISI of each individual cell. ATP concentration is 50 
pM. Black: non-oscillating cell. Circle: MDA-MB-231 cell. 
(D) When stimulated by an intermediate concentration of 
ATP (10—100 pM) the fraction of non-oscillating cells Fn 
increases with increased cancer fraction Fc at fixed total cell 
density {pr = 1200 ±200 cells/ mm^). Blue: fraction of non¬ 
oscillating NIH 3T3 cells. Red: fraction of non-oscillating 
cells including both cell types. 


(Fig. [^). Finally, we also investigate the effect of cancer 
cells on the entropy of the ISI distribution. As shown in 
Appendix B, Hjsi is higher for cells that are surrounded 
by a large number of cancer cells, and lower for cells with 
pure fibroblast neighbors. In the latter case, Hjsi also 
increases as the number of nearest neighbors decreases. 
These findings imply that reduced cell-cell communica¬ 
tion increases the entropy of the ISI values, even at the 
local level of a cell’s microenvironment, which is consis¬ 
tent with the effects seen in Fig. 0). Taken together, we 
conclude that the calcium dynamics of individual cells are 
strongly regulated by the degree of gap junction commu¬ 
nication inside the cell monolayer. 
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IV. DISCUSSION 

We have characterized the collective calcium dynamics 
of multicellular networks with varying degrees of cell-cell 
communication when they respond to extracellular ATP. 
These networks consist of cocultured fibroblast (NIH 
3T3) and cancer (MDA-MB-231) cells. We have revealed 
their properties using a combination of controlled exper¬ 
iments and stochastic modeling. We have found that 
increasing the ATP stimulus increases the propensity for 
cells to exhibit calcium oscillations, which is expected at 
the single cell level. However, we have also found that 
increasing the cell density alone, while keeping the stimu¬ 
lus fixed, has a similar effect, revealing a purely collective 
component to the sensory response. Modeling suggests 
that this effect is due to an increased degree of molecular 
communication between cells. In line with this predic¬ 
tion, we have found that increasing the fraction of cancer 
cells in the monolayer reduces the oscillation propensity, 
as cancer cells act as defects in the communication net¬ 
work. 

Typical physiological concentrations of extracellular 
ATP are tens of nM and larger [34ll36] . but hundreds of 
/iM has been associated with disease This suggests 
that healthy tissues in vivo are restricted to the regimes 
of Fig. and B, and not C (200 /iM). In the regime 
of Fig. 1^, the calcium dynamics encode the stimulus 
strength. In the regime of Fig. [^, the calcium dynam¬ 
ics encode two pieces of information simultaneously: the 
stimulus strength and the cell density. One may won¬ 
der why a system would evolve to encode two pieces 
of information in the same quantity. It is tempting to 
speculate that there is a benefit to having the sensory 
response of the system be a function of the population 
density, as in quorum sensing [4] or so-called dynamical 
quorum sensing based on collective oscillations |38ti^ . 
The particular benefit of such a strategy for these cells 
is unclear. In a similar vein, one may wonder whether 
it is possible for cells to deconvolve the two pieces of in¬ 
formation from the calcium signal using a downstream 
signaling network. Such “multiplexing” has been shown 
to be possible with simple biochemical networks lU, al- 
though the ways in which dynamic information is stored 
in, and extracted from, cellular signals is a topic of on¬ 
going research nans]. 

Our results suggest that the dependence of the calcium 
response on both sensory and collective parameters per¬ 
sists despite significant cell-to-cell variability. Indeed, we 
have found that certain measures are robust to variabil¬ 
ity, such as the oscillation propensity and the entropy of 
the ISI distribution, while others are not, such as spatial 
correlations in the ISI and its dependence on the ATP in¬ 
put (frequency encoding). This implies that traditional 
measures of information encoding in calcium dynamics, 
such as frequency encoding m uni ns m US], may have 
to be rethought in contexts where cell-to-cell variability is 
pronounced. It is becoming increasingly understood that 
variability is common in cell populations, and recent ex¬ 


amples suggest it may even be beneficial. For example, 
recent studies in a related system (NF-/^B oscillations in 
fibroblast populations) also found a large degree of cell- 
to-cell variability m and demonstrated that this vari¬ 
ability allows entrainment of the population to a wider 
range of inputs [45] . 

In our model, the transition from the non-oscillatory to 
the oscillatory regime occurs due to a saddle-node bifur¬ 
cation, a critical point in parameter space where the num¬ 
ber of dynamical fixed points changes (see Appendix C). 
This transition is broadened by intrinsic noise and cell-to- 
cell variability into a critical “region”, while cell-cell com¬ 
munication causes the oscillation propensity to depend 
on cell density within this region (Fig. |3}\). Nonetheless, 
the underlying mechanism remains the critical dynamical 
transition. Our finding that this region is broad, and our 
suggestion that it may be of some functional use for the 
system, resonates with recent studies that have argued 
that biological systems are poised near critical points in 
their parameter space [461148] . However, the connection 
between dynamical criticality, as in our model, and crit¬ 
icality in many-body statistical systems, remains to be 
fully explored. 

Gap junctional communications exist among many 
types of cells. Therefore, our results may have far- 
reaching implications for other biological model systems 
such as neuronal networks or cardiovascular systems. It 
will be interesting to explore the extent to which distinc¬ 
tions in the calcium dynamics in these systems originate 
from differences in their degrees of cell-cell communica¬ 
tion. 


V. MATERIALS AND METHODS 

A. Fabrication of Microfluidic Device 

See Appendix A for detailed device fabrication and 
characteristics. 


B. Cell Culture and Sample Preparation 

NIH 3T3 and MDA-MB-231 cells were cultured in 
standard growth mediums (Dulbeccos modified Eagle 
medium (DMEM) supplemented with 10% bovine calf 
serum andl% penicillin and DMEM supplemented with 
10% fetal bovine serum, 1% penicillin, and 1% non- 
essential amino acids respectively). To prepare samples, 
cells were detached from culture dishes using TrypLE Se¬ 
lect (Life Technologies) and suspended in growth medi¬ 
ums before pipetted into the microfluidics devices and 
allowed to form monolayers. If MDA-MB-231 cells were 
the dominant species (a fraction greater than 50% of all 
cells), they were first allowed to attach the glass bottom 
of the microfluidics devices. Red fluorescent tag (Cell- 
Tracker, Life Technologies) was then applied and subse¬ 
quently washed with growth medium. Einally, NIH 3T3 


cells were injected into the device so that the desired cell 
density (^1000 cells/mm^) was reached. If NIH 3T3 cells 
were the dominant species, they were allowed to attach 
the glass bottom of the microfluidics devices first. MDA- 
MB-231 cells already loaded with CellTracker were then 
injected into the devices to reach the desired cell density. 
After incubating the microfluidics devices containing cell 
monolayers overnight, fluorescent calcium indicator was 
applied (Fluo4, Life Technologies) making the samples 
ready for imaging. 


suing calcium dynamics of cells of various densities on a 
3-by-3 grid (and 7-by-7 grid for Fig. using the Gille¬ 
spie algorithm |29] with tan-leaping 49]-[52]). We include 
cell-to-cell variability by sampling all model parameters 
from distributions that are uniform in log space, which 
varies fold-change up to a maximal value. For further 
details, see Appendix C. 


C. Fluorescence Imaging and Image Analysis 

Fluorescence was detected using a inverted microscope 
(Leica DMI 6000B) coupled with a Hamamatsu Flash 2.8 
camera. Movies were taken at a frame rate of 4 frame/sec 
with a 20x oil immersion objective. Image analysis and 
data processing were performed in MATLAB. (Details in 
SI). 

D. Stochastic model of collective calcium dynamics 

We extend the model in [Ml El to include intrinsic 
noise and cell-cell communication. We simulate the en- 
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Appendix A: Additional information of experimental setup 

1. Fabrication and characteristics of microfluidic device 

The organic elastomer polydimethylsiloxane (PDMS, Sylgard 184, Dow-Gorning) used to create the microfluidic 
devices is comprised of a two part mixture - a base and curing agent - that is mixed in a 10:1 ratio, degassed, and 
poured over a stainless steel mold before curing at 65°G overnight. Once cured, the microfluidic devices are cut from 
the mold, inlet/outlet holes are punched, and the device is affixed to a No. 1.5 coverslip via corona treatment. Figures 
[^-G provide schematics of the design of the device as well as inner dimensions of the flow chamber. 

To characterize the stimuli profile in the flow chamber, we use fluorescein (Sigma Aldrich) to replace ATP as our 
probe and record the fluorescent intensity in the same flow condition as in the ATP sensing experiments. Figure 
shows the normalized fluorescent intensity as a function of time from two devices. We estimate the time takes to fill 
the field of view with chemical stimuli is about 1 second. 


2. Obtaining individual cell calcium dynamics 

Application of GellTracker Red GMPTX (LifeTechnologies) to breast cancer cells and the calcium binding dye 
Alexa Fluo4 (LifeTechnologies) to both breast cancer cells and fibroblast cells within the microfluidic device allows 
for the differentiation of the two cell types. To do so, successive images are taken of the two fluorophores and first 
analyzed in Imaged (imagej.nih.gov). Overlaying the two images generates a composite image that determines the 
network architecture, as is shown in Figure lA of the main text. The centroids for each individual cell are manually 
determined and a dot is placed at the center of each cell as shown in Fig.[^ Using the dotted image as a mask, intensity 
is averaged over the ^40 pixels comprising the dot for each cell and for all acquired frames in the experiment. 


3. Quantification of gap junctional intercellular communication 

In order to quantify the gap junction communications for the two cell types we have used (NIH 3T3 and MDA-MB- 
231), we perform FRAP (fluorescent recovery after photobleaching) experiments to measure the effective diffusion 
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FIG. 6: Design of the microfluidic device and stimuli arrival time (A) Top-view of microfluidic device showing glass slide, 
PDMS device, inlet/outlet, and connecting microfluidic chamber. (B) Side-view of microfluidic device indicating flow direction 
in red. (C) Cross-sectional view of device providing dimensions of flow layer and attached fibroblast (Green) and breast cancer 
(Red) cell monolayer. (D) Temporal prohles of chemical stimuli in the flow chamber evaluated using fluorescein. Two devices 
are tested (red and blue). Vertical lines correspond to arrival time and the times when half-maximum intensity reached. 


coefficients in a cell monolayer for each cell type, (some experimental details, eg dye used frame rates etc). As 

shown in Fig.|^, at each time point t, we represent the average concentration of fluorescent dye in the bleached area 

as Ci{t) and non-bleached area as Co{t). The concentration of the dye is proportional to the recorded intensity as 

C{t) = al{t), where a is a constant. We approximate the time evolution of Co{t) and Ci{t) based on the following 

assumptions: (a) there is an overall bleaching rate P while imaging, affecting both Co{t) and Ci{t). (b) Ci{t) get an 

influx due to a concentration gradient approximately proportionally to -Pe// where P^e// is the effective 

2 

diffusion coefficient, (c) the cell monolayer has a uniform thickness of h. We have 

ACo{t) = Co{t + At) - Co{t) « -pCo{t)At (Al) 

ACi{t) = Ci{t + At) - Ci{t) 

^ ^Deff[Co{t+^)-Ci{t + ^)]\At-/3Ci{t)At (A2) 

The above equations can be simplified and we have 

p AJ.W - 

In our experiment, we take image at every 30 secs, therefore At is set to be 1 minute. The recorded time series 
/o(t), Io{t) from all experiments are plugged into the above equation to obtain a distribution of the effective diffusion 
coefficient D^ff as shown in Fig. [^. 
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FIG. 7: Determination of single cell centroids. (A) Composite image showing all cells (NIH 3T3 and MDA-MB-231, green), 
MDA-MB-231 (red), and manually determined centroids of NIH 3T3 cells (blue). (B) Composite image showing all cells (NIH 
3T3 and MDA-MB-231, green), MDA-MB-231 (red), and manually determined centroids of all cells (NIH 3T3 and MDA-MB- 
231, blue). 


Appendix B: Additional information of Inter-Spike Interval (ISI) analysis 


1. Inter-Spike Interval Determination 


As calcium spike trains are unsuitable for Fourier analysis due to their varying amplitudes and irregular periods, the 
inter-spike interval (ISI) has been used as a suitable substitute for determining how the period of calcium oscillation 
affects information encoding within and amongst cells in the network. To determine the ISFs from a single cell’s 
calcium dynamics i7(t), we first smooth the time series R{t) to remove outlying data points [53], to obtain Rs{t). We 
then find all local peaks {ti} in Rs{t) using the following criteria: a local peak must be higher than its neighboring 
valleys (local minimums) on both directions by a certain value, which we refer to as the sensitivity parameter; the 
first and last frames are excluded as peaks or valleys. Third, we obtain ISFs {6ti} from the distances of peaks 
{6ti = — ti\}. Finally, to obtain the true ISFs corresponding to calcium oscillations, we only keep {6ti} fall in 

the range of 5 seconds and 150 seconds. This final step excludes about 5% of ISFs. 

Figure and[^ provide two examples of the peak finding algorithm’s ability to determine peak locations for a 
fibroblast cell calcium profile (A) and cancer cell calcium profile (B). The sensitivity parameter used for all experiments 
and in Figure]^ andj^ is | the standard deviation of the normalized intensity aj as determined empirically. Figure 
and[^ show the algorithms diminished ability to determine peaks and to detect spurious peaks when the sensitivity 
parameter is set to aj and ^ respectively. 
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FIG. 8: Estimation of the effective diffusion coefficients mediated by gap junctional communication in a cell monolayer. (A) 
Cell monolayer loaded with fluorescent dye (CellTracker Red CMPTX, LifeTechnologies) has been partially bleached. (B) 
The FRAP experiments can be quantified by the average intensity of the bleached area (A) and non-bleached area (/o). (C) 
Normalized intensity prohles of a typical experiment. Ipb is the intensity of the bleached area while being bleached. A and lo 
are intensity prohles during recovery. (D) Histogram (normalized to probability) of the diffusion coefficients for NIH-3T3 cells 
and MDA-MB-231 cells. 


2. Spatial map of ISI in multicellular network 

We have shown a map of ISI in Fig. 5C of the main text. To generate such a map, we first calculate the ISFs 
corresponding to each individual cell as shown in the above. We then calculate the average ISI of each cell. Fig. 
[Tq| shows two more examples of the ISI spatial map. Figure pT| A corresponding to a map from a experiment with 
20/iM ATP concentration and Fig. corresponding to a 50/iM ATP concentration. Black regions within the maps 
indicate cells that did not oscillate in the given experiment and a red circle indicates the location of a cancer cell. It 
is apparent from the figure that cancer cells are more likely to be non-oscillatory than fibroblast cells. This is also 
shown in Fig. El which shows ISI distributions of both cell types for a typical experiment. 


3. Entropy of ISI distributions 

Interestingly, the network architecture also regulates the spectrum of ISI. This can be demonstrated from the 
perspective of each individual cell’s microenvironment. We have grouped NIH 3T3 cells based on their number of 
nearest neighbors (A^nn)? ^.s well as the fraction of cancer cells in their N^n neighbors. For each group at a fixed 
ATP concentration (50 jiiM), we have computed the differential entropy of interspike intervals Hjsi. The differential 
entropy is computed using the k-nearest neighbor method, as described in Appendix D. As shown in Fig. Hjsi 
is higher for cells that are surrounded by a large number of cancer cells, and lower for cells with pure fibroblast 
neighbors. In the latter case, Hjsi also increases as the number of nearest neighbors decreases. These results show 
that the spectrum of ISI is regulated by cell-cell communications both locally in the immediate microenvironment of 
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FIG. 9: Intensity peak determination for single NIH 3T3 and MDA-MB-231 cell profiles. (A) Labeled peaks for a NIH 3T3 
cell. Sensitivity parameter set to | of the standard deviation of the time series. (B) Labeled peaks for a MDA-MB-231 cell. 
Sensitivity parameter set to | of the standard deviation of the time series. (C) and (D) Same MDA-MB-231 intensity profile 
in (B) but with sensitivity parameter set to be the standard deviation and | of the standard deviation of the time series 
respectively. 


a cell and globally throughout the whole defective multicellular network. 


4. Statistics of FV, the fraction of non-oscillating cells 

As shown in Figs. 4 and 5 of the main text, the fraction of non-oscillating cells systematically depends on 
the network architecture. Fig. shows pairwise statistical comparison of the results. The results are obtained by 
ANOVA followed with Fisher’s least significant difference procedure, implemented in MATLAB. 


Appendix C: Theoretical Model 
1. Deterministic Model 

We model the dynamics of calcium release following Tang and Othmer 120 ]. Their model is deterministic and 
valid for a single cell. First we describe the behavior of their model. Then we generalize it to include stochastic 
effects. Next, we generalize it further to include cell-cell communication. Finally, we add cell-to-cell variability. 
This procedure allows us to systematically explore the effects of stochasticity, communication, and variability on the 
calcium dynamics. 
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FIG. 10: Spatial map of cell-average ISI of randomly selected experiment with (A) 20/iM ATP. (B) 50/iM ATP. 


ATP triggers the production of IP 3 within a cell. IP 3 binds to receptors on the endoplasmic reticulum (ER). The 
ER is an internal store for calcium ions. When the receptors are active, they act as channels that allow passage 
of calcium between the ER and the cytoplasm of the cell. Receptor activation requires binding of both IP 3 and 
cytoplasmic calcium, such that calcium triggers its own release. However, an excess of calcium deactivates receptors, 
closing the channels. These positive and negative feedbacks, respectively, lead to oscillations and other rich dynamics. 

The model of Tang and Othmer is based upon the following reactions. 


90 ^ 

c, 

Cs + R4^C + i?4, 

c^c,, 

(Cl) 

^90 


UJi 


kil 
_ t _ k 

- 

k-i 

R3, 

R3 + C \ --R4, 

k-2 

R4 T C ^ - -R5, 

k-3 

(C2) 


Eqn. |C1| contains the reactions describing calcium transport between the ER and the cyto plasm , while Eqn. |C 2 | 
contains the reactions describing the sequential (un)binding of molecules to receptors. In Eqn. Cl, the first reaction 
describes leakage between the calcium in the ER (Cs) and calcium in the cytoplasm (C); v is the ratio of the ER 
volume Vs to the cytoplasm volume V. The second reaction describes transport of calcium through the receptor 
channels: receptors with both IPS and calcium bound (R 4 ) are active channels. The last reaction describes the 
pumping of calcium from the cytoplasm to the ER; / denotes the nonlinear pumping propensity. In Eqn. |C2| the 
first reaction describes the binding of IPS to bare receptors {R 2 ), creating the complex R 3 ; / is the concentration of 
IPS, which is treated as a parameter in their model (we later generalize this feature, making I a dynamic variable 
and its production rate a the control parameter, since this is more consistent with the experimental setup). The 
second reaction describes the subsequent binding of cytoplasmic calcium to the complex, which activates the complex 
(R 4 ). The last reaction describes the further binding of cytoplasmic calcium to another site on the complex, which 
deactivates the complex (R 5 ). Only receptors in the R 4 state serve as active channels. 
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FIG. 11: ISI characteristics of a typical experiment. ISI event counts normalized by number of cells. Blue: statistics of only 
NIH 3T3 cells, Red: statistics of only MDA-MB-231 cells. 
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The deterministic model describes the dynamics of the mean numbers of molecules. The dynamics follow directly 
from Eqns. |C1| and |C2[ Introducing n and Ug as the numbers of calcium molecules in the cytoplasm and ER, 
respectively, and rrii as the numbers of receptors in state i G {2, 3,4, 5}, the means (denoted with an overbar) obey 


dn 

dt 

dfig 
dt 
dm 2 
dt 

dfhs 

dt 

dfri 4 

dt 

dffi^ 

dt 


goUs - ygon - z/tti 
Mfh^n + k- 2 ^A - 


—2 + 71^5^4 - v^\nm^ 
+ 7r2 


-goUg + lygon + utti — 2 “ + z^7i^^4, 

+ 7r2 


-kilfh2 + k-iffis, 
kilfh2 — k-ifhs — \2Tri^n + /c_2m4, 
\2fh^n — k-2kfi4^ — Xsffi^n + k-^ffi^^ 
X^ffi/^n — k-^ffi^. 


(C3) 

(C4) 

(C5) 

(C6) 

(C7) 

(C8) 


The pumping propensity is a Hill function with coefficient 2 and parameters tti and 7r2. 

Two quantities are conserved in the model: the total number of calcium molecules no and the total number of 
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FIG. 12: The spectrum of inter-spike intervals (ISI) as characterized by Hisi is regulated by the microenvironment of cells. For 
hxed ATP concentration (50 /xM), calcium dynamics and the associated ISI of individual hbroblast cells are grouped based on 
the microenvironment of the cells. The microenvironment is characterized by the number of nearest neighbors (A^nn), among 
which a fraction of Fnnc are cancer cells. The color map represents Hi si (in bits) at varying Nnn and Fhnc, obtained by 
interpolating scattered data points with a Gaussian kernel. 
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FIG. 13: Statistical analysis of Fn with varying cancer fractions (A) and cell densities (B). For both cases, the middle two 
groups do not show statistically signihcant differences. Therefore we have treated the middle two groups as a single group and 
performed ANOVA analysis. Notice that the blue bars, which represent non-oscillating hbroblast cells are not analyzed for p 
values. 
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Deterministic parameters 

Stochastic parameters 

jy = 0.185 

Co = 1.56 /jM 
go = 0.025 s“^ 
gi = 36 s“^ 
pi = 54 

P2 = 0.03 jjM 

ki = 120 

k2 = 150 
ks = 0.18 
k-i = 96 s“^ 
k-2 = 18 s“^ 
k-3 = 0.018 s-^ 

V = 200 

mo = 200 

ly = 0.185 

no = (1 + i^)coV = 2.2 X lO'^ 
go = 0.025 

71 = gi/mo = 0.18 

TTi = piV = 6.5 X 10® 

— P 2 V — 3600 

Ai = ki/V = 0.001 s“^ 

A2 = fe/E = 0.0013 s-i 

A3 = fc3/E= 1.5 X 10“® s-i 

fc-i = 96 

k-2 = 18 s“^ 

k-3 = 0.018 s-i 

a (varied) 

y = 0.058 s“^ 

h = 0.025 s“^ 


TABLE I: Parameters and their values. Deterministic values are taken from Tang and Othmer [20]. Stochastic values follow 
from the deterministic values, given additional choices for the cytoplasmic volume V and the total number of receptors mo- 
These two parameters, as well as the additional parameters a, /li, and h, are estimated as described in the text. 


receptors mo, 


no = n + ns+m 4 + 2 m 5 , (C9) 

mo = m 2 +m 3 +m 4 + m 5 . (CIO) 


The conservation relations reduce the number of dynamical equations from six to four. Specifically, we eliminate n^ 
and m 4 . Furthermore, we assume that the number of calcium molecules is much larger than the number of receptors, 
no ^ mo- In this limit, introducing c = (n + m 4 + 2fh^)lV as the concentration of calcium in the cytopla sm (both 
unbound and bound to receptors) and Xi = mijm^ as the fraction of receptors in each state, Eqns. C3][C8 become 


dc 

dt 

dX 2 

dt 

dxs 

dt 

dx^ 

dt 


(1 + p){go + giX4){co - c) - i/pi 



-kiIX2 + k-iX3, 


kiIX2 - k-iX3 - k2X3C + k-2XA, 


k3XAC - k-3X5, 


(Cll) 

(C 12 ) 

(C13) 

(C14) 


where 2:4 = 1 - a ;2 - 2:3 - 2 : 5 , gi = 7 imo, pi = ni/V, p 2 = tt 2 /V, k 2 = A 2 y, fca = M V, and cq = no/(V + K) is 
the average calcium concentration in both the cytoplasm and the ER. Eqns. |C11||C14 reproduce Eqn. 4 in |20] . The 
parameter values used in 120], which we also use here, are given in Table | 

Tang and Othmer recognize that the (un)binding of calcium to the deactivating site is slow compared to the other 
(un)binding reactions, ks {ki,k 2 } and ks <C {fe-i, k - 2 }. This permits the quasi-steady-state approximations 
dx 2 ldt ~ 0 and dx^/dt ~ 0, which reduce Eqns. C11|C14 to two equations. In dimensionless form they read 


dnr 

, = a 3 {l-x) + a 2 {l-x) 

dr 


x{l - y) 


x + j3i[l + j5o{I)] x'^+al 


I32x‘^{ 1 - y) 

X + Pi[\. + Poil)] ’ 


(CIS) 

(C16) 


where x = c/cq, y = x^, t = k_3t, e = k_3Co/vpi, ai = (1 + vyogolvpi, a2 = {I + = P2/C0, 

/3o(/) = k-i/kil, pi = k- 2 /k 2 Co, and /?2 = k-sco/ks. 
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FIG. 14: The deterministic model (Eqns. C15 and C16) exhibits four dynamic regimes: (A) monostable, with low cytoplasmic 
calcium, (B) excitable, (C) oscillatory, or (D) monostable, with high cytoplasmic calcium. In the last regime (D), near the 
boundary, oscillations can also be supported for certain initial conditions. (E) These regimes are accessible by tuning the IP3 
concentration /. The letters A-D in panel E correspond to the dynamics seen in panels A-D. The critical values /i, , and 

separating the regimes (dashed lines) are determined by changes in the number or stability of fixed points, which follow from 
a standard linear stability analysis. Parameters are listed in Table ^ 


Eqns. C15 and C16 exhibit four dynamic regimes. As shown in Fig. 14 


the system can be (A) monostable, 
with low cytoplasmic calcium, (B) excitable, (C) oscillatory, or (D) monostable, with high cytoplasmic calcium. In 
the last regime (D), near the boundary, oscillations can also be supported for certain initial conditions. All four 


regimes are accessible by tuning the IP 3 concentration I (Fig.[l4p). Transitions between regimes occur at the critical 
concentrations , and /|. Of particular interest is the transition from the excitable to the oscillatory regime (/I )5 

since these are the dynamics exhibited by cells in the experiments: transient pulsing or sustained oscillations. 


Stochastic Model 


Stochastic dynamics are simulated using an adaptive tan-leaping method [49115^ . which is a computationally- 
efficient approximation of the Gillespie algorithm [29]. The simulation accounts for the facts that molecule numbers 
are integer-valued and that reactions fire at random, exponentially-distributed times. The reaction propensities follow 
directly from Eqns. C 1 |C 2 , and have the same forms as the individual terms in Eqns. C3||C8| (without the overbars). 
The stochastic model parameters are related to the deterministic model parameters by the cytoplasmic volume V 
and the total number of receptors mo; see Table [I| Together these factors determine the noise in the system via the 
molecule numbers tiq = {1 u)cqV and mo- Since relative fluctuations decrease with molecule number, large no 
corresponds to low noise, while small no corresponds to high noise. Typical trajectories are shown in Fig. 

We estimate the cell volume from the experiments. An upper bound is obtained from the measurements of the cell 
density p, since the cell area cannot be larger than 1/p. The maximal density p = 2500 mm“^ gives the sharpest upper 
bound, for an area of 400 pm^, or a length scale of 20 pm. Since this is an upper bound, we assume a length scale 
range of 10—20 pm. For computational efficiency, we perform simulations assuming 10 pm, and taking a monolayer 
width of ^2 pm, this leads to a cell volume of ^ = 200 pm^. Absent an experimental estimate for the number of 
receptors mo, we take mo = 200 . 

In the stochastic model, we no longer treat I (the concentration of IPS) as a parameter; instead, we treat the 
number of IPS molecules as a stochastic variable. The control parameter is then the production rate of IPS, a, which 
is assumed to be proportional to the ATP concentration in the experiments. IPS is degraded at a rate p, meaning 
that the stochastic reactions simulated are Eqns. icn^ with the first reaction of Eqn. |C 2 | replaced with 


1^0, 


i?2 + I ^ -^ i?3- 


k-i 


(C17) 


The degradation rate of IPS has been measured to be p = 0.058 s' -1 m, and we vary a in the range 


-lO^-lO"^ s-i 
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FIG. 15: Noise expands the oscillatory regime. (A) Color-co ded bar graph showing the dynamica l regi mes as a function of I 
determined by the reduced deterministic model (Eqns. C3|C8 top), the deterministic model (Eqns. C15 and C16) for a chosen 
set of initial conditions (middle), and the stochastic model using the same initial conditions (bottom). (B-D) Typical dynamics 
(top) and scaled power spectra (bottom) for the deterministic model (red) and stochastic model (black) are shown for the 
/-values denoted by the dashed lines labeled B-D in A. 


such that the mean IPS level is (a//i ^ 10^ — 10^ molecules per cell, or 0.1 — 1 /iM. 


3. Cell-to-Cell Communication 

We introduce communication into the system by allowing calcium molecules to diffuse between neighboring cells, 
simulating communication via gap junctions in the experiments. Doing so introduces new reactions into our model: 

Cij ^ C'i+i.j, Cij ^ C'i.j+i, (C18) 

h h 

where i,jf index cell lattice sites, and h is the hopping rate of calcium ions from cell to cell. We estimate h from 
the experimentally measured diffusion coefficient D ^ 1 /rm^/s (see Fig. ID of the main text). Since in our model, 
ions are hopping between cells on a two-dimensional lattice, h~^ corresponds to the expected time ^ j^D to diffuse a 
cell-to-cell distance Using ^ ^ 10—20 jam as above, we find h ^ 0.01—0.04 s“^, and therefore we use h = 0.025 s“^. 

Our calculations are done on a square lattice, where each lattice site is either empty or contains one cell, and 
therefore each cell can have up to four neighbors with which to communicate. Density is varied by changing the 
number of empty lattice sites. For each chosen density, we sample over individual realizations, in which cell locations 
are assigned randomly. Thus, the statistics encompass many possible spatial distributions of cells for each density. 
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4. Effects of noise and communication 


Once we add noise, we can no longer distinguish the dynamic regimes (monostable, excitable, oscillatory) using the 
deterministic criterion, namely the number and stability of the fixed points (Fig. [T^). Instead, we define a criterion 
based on the trajectories themselves, similar to the procedure for the experimental trajectories. We simulate the 
system for a given time T, where T = 400 s as in the experiments. The criterion is then: 

• If the calcium level spikes zero times in T, the dynamics are monostable (non-oscillatory), 

• If the calcium level spikes one time in T, the dynamics are excitable (non-oscillatory), and 


• If the calcium level spikes two or more times in T, the dynamics are oscillatory. 

When analyzing trajectories to determine ISI values, excitations separated by less than t = 5 s are excluded, as are 
excitations separated by more than t = 200 s, just as in the experimental analysis. We validate the above criterion 
using a complementary approach: we also compute the power spectrum of the trajectories. For periodic trajectories, 
the power spectrum exhibits a peak at the oscillation frequency. As shown in Fig.[T^-D, trajectories that are periodic, 
and thus pass the above criterion, indeed given peaked power spectra as expected. 

We find that noise broadens the boundaries between dynamical regimes (Fig. 2B of the main text), and that noise 
shifts the boundaries between regimes (Fig. 15 4, compare the bars in the middle and bottom rows). In particular, 
the oscillatory regime is expanded at both ends. At low values of the control parameter, noise promotes oscillations 
where the deterministic dynamics are excitable. At high values of the control parameter, noise promotes oscillations 
where the deterministic dynamics are mono-stable. Both effects have been observed in stochastic models of similar 
excitable systems, and are likely due to the ability of noise to cause repeated excitations and prevent damping of 
oscillations, respectively [55] . 


5. Cell-to-Cell Variability 

The model predicts frequency encoding, but experimental measurements show no trend in ISI values as a function 
of ATP concentration (Fig. IE). We thus introduce cell-to-cell variability into the model as a possible explanation for 
the lack of observed frequency encoding. 

Variability is introduced into our model by allowing parameters for each cell to be drawn from a distribution. All pa¬ 
rameters {mo, 'Uo, ^0 5 tti, 7 r 2 , Ai, A 2 , As, /c_i, /c_ 2 , A^-s, /i, h} are allowed to vary from their baseline values listed 

in Table|^ Parameters are sampled from uniform distributions in log space such that logio(x) G [logio(x/F),logio(xF)], 
where x is the mean value for the parameter x, and F is a parameter that changes the strength of variability. Therefore 
each parameter is sampled from a distribution such that x G [x/F^ FF \. 

Fig. [T6| shows the ratio of the standard deviation in mean ISI values over their mean as a function of F. As F is 
increased, this ratio increases until saturation occurs at F ^ 2 . Also shown is box plot of the experimental values 
of this ratio, which was obtained from twelve experiments with similar density. By the saturating value F = 2, the 
theoretical values he within the experimental range; therefore we use the value F = 2 to obtain the results that include 
variability in the main text. 


Appendix D: Nearest neighbor cross-correlation analysis and computation of differential entropy 

1. Definition 

We have previously demonstrated the effectiveness of using nearest neighbor cross-correlation functions (Cnn) to 
characterize the calcium dynamics of collective chemosensing m- As a result, we first study how network architecture 
may affect the spectrum of Cnn- For any cell i, we represent its calcium dynamics (fluorescent intensity profile) as 
Fi(t), and define Cnn as: 

CNNij = C{t = 0)ij \d..=i ■ (Dl) 

To evaluate we numerically differentiate the response curve Ri{t) using the five-point stencil method; 
is the time average. Note that is the variance of Ri{t)^ which normalizes C{r)ij to be dimensionless. Based on 
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FIG. 16: The effect of changing the strength of variablity parameter, F. The ratio of the standard deviation to the mean of 
the ISI, as a function of F. The ratio saturates above F ^ 2. The range of values for experiments with similar density is shown 
at right. 


TABLE II: Statistics of nearest neighbor pairs 


[ATP] 

10 pM 

20 pM 

50 pM 

100 pM 

Total 

# of cells 

7,055 

4,245 

13,078 

5,062 

30000 

# of nearest neighbor pairs 

40,596 

24,462 

75,092 

32,222 

172,372 


Delauney triangulation of the multicellular network, topological distance D can be defined for each cell pair where 
D = 1 for nearest neighbors. The mean nearest neighbor cross-correlation function Cnn is obtained by averaging 
C^NNij over all nearest neighbor pairs i^j. 


2. Statistics of nearest neighbor pairs 

By using Delauney triangulation coupled with the manually determined cell centers (see Fig.for example), we are 
able to determine the connectedness of the monolayer network. In particular, we focus on the communication between 
nearest neighbor cells to reflect proper intercellular communication which correspond to a topological distance of 
D = 1 in the triangulation. The total number of nearest neighbor pairs and cells analyzed for each ATP concentration 
are given in table |n[ 


3. Hcnn of collective calcium dynamics 

We compute the spectrum of CnNij using differential entropy, a scalar value that characterizes the randomness of 
a set of observables. The differential entropy is defined as H = — f p(x)p(x)dx, where p(x) is the probability density 
of a continuous variable x. Higher values of differential entropy are associated with broader distributions; this allows 
us to directly compare the spectrum of nearest neighbor cross-correlations. To avoid biasing the entropy calculation 
due to binning of data, the differential entropy is determined by using the /c-nearest neighbor algorithm [SHIISI] , as 
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FIG. 17: Differential entropy of nearest neighbor cross-correlations (Hcnn) as function of cell density at fixed fraction of cancer 
cells (15%±4%) 


described below. As shown in Fig. increasing cell density decreases the differential entropy of Cnn^ suggesting 
that higher levels of network connectivity will suppress the fluctuations in the nearest neighbor cross-correlations. 


4. Evaluations of differential entropies 

In order to characterize the spectrum of nearest neighbor cross-correlations and inter spike intervals of the collective 
calcium dynamics, we have implemented a k-nearest neighbor method to calculate the differential entropy of a set of 
random variables. We first make benchmark tests using random variables drawn from Gaussian distributions. Since 
most experimental datasets contain more than 1500 elements, we compare the differential entropy calculated from 
1500 Gaussian random variables with respect to the analytic results. As shown in Fig. [T^, numerical calculations 
agree with exact values very well. The performance of k-nearest neighbor method depends on the sample size. To 
evaluate the dependence, we employ k-nearest neighbor method on random variables drawn from standard normal 
distribution with varying sample sizes. As shown in Fig. [T^, the relative error is around 7% for sample size equals 
to 100, and quickly decay to less than 1% when sample size is greater than 500. 

In calculating the differential entropy H^iff using k-nearest neighbor method, every element in a dataset contribute 
to the final result. In order to estimate how much Hdiff fluctuates as a result of resampling, we would like to calculate 
Hdiff from many datasets each of which contains distinct elements drawn from the same probability distribution. 
However, the distribution of experiment observables are unknown. In order to perform resampling, we employ statis¬ 
tical bootstrap [58]. As a particular example, we take a randomly selected set of nearest neighbor cross-correlations 
and calculated the differential entropy of 1000 bootstrap resamplings. Fig. [T8|C shows the probability distribution of 
the results. As can be seen, the fluctuation is small, and is typically less than 2%. 
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FIG. 18: Evaluation of differential entropy using k-nearest neighbor method. A: Differential entropy Hdiff for normal distri¬ 
butions of varying standard deviations a and zero mean. For each cr, 1500 random variables are drawn and k-nearest neighbor 
method is used to calculate the corresponding differential entropy. 100 realizations (redraw the random variables) are evaluated 
for each a to calculate the mean and standard deviation of Hdiff shown as the error bars. Dashed line is the analytic result. 
B: Relative errors of differential entropy Hdiff calculated with varying sample sizes. For each data set of sample size n, n 
random variables are drawn from Gaussian distribution with zero mean and standard deviation equals to 1. 100 realizations 

are made for each sample size to calculate the relative error of Hdiff. Relative error is defined as where Hi is 

the differential entropy calculated with k-nearest neighbor method at realization z, and Ho = Ln{27re)log2e. G: Histogram of 
Hcnn from bootstraping a randomly selected experimental data set. Hcnn of 1000 bootstrap realizations are calculated and 
the corresponding histogram is normalized to probability. Dashed line indicates the mean value. 
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